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We numerically compute the helicity modulus of the three-dimensional gauge glass by Monte Carlo 
simulations. Because the average free energy is independent of a twist angle, it is expected that 
the average helicity modulus, directly related to the superfiuid density, vanishes when simulations 
are performed with periodic boundary conditions. This is not necessarily the case for the typical 
(median) value, which is nonzero, because the distribution of the helicity modulus among different 
disorder realizations is very asymmetric. We show that the data for the helicity modulus distribution 
are well described by a generalized extreme-value distribution with a nonzero location parameter 
(most probable value). A finite-size scaling analysis of the location parameter yields a critical 
temperature and critical exponents in agreement with previous results. This suggests that the 
location parameter is a good observable. There have been conflicting claims as to whether the 
superfiuid density vanishes in the vortex glass phase, with Fisher et al. [Phys. Rev. B 43, 130 
(1991)] arguing that it is finite and Korshunov [Phys. Rev. B 63, 174514 (2001)] predicting that it 
is zero. Because the gauge glass is commonly used to describe the vortex glass in high-temperature 
superconductors, we discuss this issue in light of our results on the gauge glass. 



PACS numbers: 75.50.Lk, 75.40.Mg, 05.50.+q 

I. INTRODUCTION 

The gauge glass is often used to describe the ex- 
tremely disordered vortex glass phase in granular high- 
temperature superconductors Fi^ii In addition, it is the 
simplest possible model that possesses the correct order 
parameter symmetry [XY spins with a U(l) symmetry] 
and therefore considerable numerical, as well as theoret- 
ical, work has been based on it. Nevertheless, the model 
has some limitations: it is isotropic, whereas finite mag- 
netic fields in the vortex glass phase in high-temperature 
superconductors introduce a uniaxial field anisotropy. 5 
Furthermore, the gauge glass ignores transverse screen- 
ing, an essential ingredient in superconductivity. In fact, 
the inclusion of screening destroys the gauge glass phase 
and the critical temperature drops to zero£ In addi- 
tion, the disorder is introduced via random gauge fields, 
whereas a more realistic model would probably introduce 
the disorder via random bondsi^ 

In three space dimensions the (unscreened) gauge glass 
undergoes a phase transition into a glassy phase at 
T c = 0.47(3) which has been estimated by a finite-size 
scaling analysis of the currents in the system (first-order 
derivative of the free energy with respect to an infinites- 
imal twist to the boundaries) »ii Another central quan- 
tity to study is the helicity modulus — directly related to 
the superfiuid density — which corresponds to the second- 
order derivative of the free energy with respect to an in- 
finitesimal twist to the boundaries. Motivated by contra- 
dicting predictions as to whether the superfiuid density 
p s is finite or not for a bulk vortex glass system ; 12 ' 13 we 
have studied in detail the helicity modulus of the three- 
dimensional gauge glass model. Because the simulations 
are performed using periodic boundary conditions in or- 
der to avoid finite-size effects due to small system sizes, 
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FIG. 1: (Color online) Possible phase diagram of a disordered 
superconductor (Ref. |4J) . The vortex lines are pinned by de- 
fects and the vortex solid found in a clean superconductor 
is transformed into a glass. For low disorder, there are no 
dislocations and one obtains a Bragg glass which exhibits a 
superconducting response, whereas for high disorder and low 
enough, temperatures a highly disordered vortex glass phase 
emerges (Ref. 0) . The vortex liquid remains liquid after the 
inclusion of disorder, albeit viscous and separated by a glass 
transition line H g (dotted line) from the vortex glass phase. 
H c i separates the Bragg glass phase from the Meissner phase, 
whereas H C 2 denotes the crossover from the superconducting 
to the normal phase (Ref. FToT) (mean-field calculation, dashed 
line). The melting lines between Bragg glass, vortex glass, 
and vortex liquid are denoted by Hm and are calculated tak- 
ing into account fluctuations. (Figure adapted from Ref. |j). 
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we expect that the average of the helicity modulus is 
zero. This (unphysical) result follows from twists forced 
into the system due to the used periodic boundary con- 
ditions. The problem can be alleviated by using either 
free boundary conditions or fluctuating-twist boundary 
conditions , 14 i 15 where location-dependent twists are ap- 
plied along the boundaries. The former have the problem 
that corrections to scaling are huge for numerically ac- 
cessible system sizes and the latter generally overestimate 
the stiffness of the system. 

Our results with periodic boundary conditions show 
that the mean of the helicity modulus does not have crit- 
ical scaling, thus preventing us from determining the lo- 
cation of the glass transition. In addition, strong fluctu- 
ations in the mean develop below the transition temper- 
ature. In contrast, the typical (median) value shows the 
correct scaling, i.e., the data cross at the known value 
of the transition temperature. Furthermore, we demon- 
strate that the data arc well fitted by extreme-value dis- 
tributions, commonly used in economics studieS ) 16 ' 17 i 18 
explaining the strong fluctuations in the mean as well 
as providing detailed estimates of the critical parameters 
via a finite-size scaling analysis of the location parame- 
ter of the helicity modulus distributions (most probable 
value). 

The gauge glass has attracted interest in the context 
of disordered high-temperature superconductors. The 
phase diagram of such disordered type-II materials ex- 
hibits numerous phases 4 -! 2 -^ (see Fig. [T]) and, despite 
intensive research, several features of the different phase 
diagrams of high-temperature superconductors in a mag- 
netic field remain to be fully understood. In the presence 
of disorder the Abrikosov lattice ceases to be perfect. 
For low enough temperatures and fields (above H c \), a 
Bragg glass emerges 4 characterized by the appearance 
of algebraic Bragg peaks when spectroscopic measure- 
ments are performed, even though the flux line lattice 
is disordered. These are lost when either the disorder or 
the temperature is increased. Increasing the temperature 
yields a vortex liquid, where flux lines fluctuate subject 
to a random disorder potential. Increasing the disorder 
yields a vortex glass phase*^ where the order of the flux 
line lattice is destroyed and the vortex lines are strongly 
pinned to the impurities with no spatial order - simi- 
lar to standard spin glasses 2 ^ (see Fig. [T]) . Signatures of 
a vortex glass phase have been observed experimentally 
by studying the I—V characteristics of YBCO samples. 
There the vanishing of the resistivity for fields larger than 
H c i has been interpreted as the onset of a vortex glass 
phasej 21 ' 22 ' 23 Furthermore, a recent microscopic study 2 ^ 
on LSCO samples using muon-spin rotation experiments 
and small-angle neutron scattering shows the transition 
from a Bragg to a vortex glass phase. 

Fisher et alw^ argue that the superfluid density p s is 
finite in the thermodynamic limit, with finite-frequency 
corrections that only vanish as an inverse power of ln(w). 
In contrast, KorshunoA*^ recently predicted that the su- 
perfluid density tends to zero logarithmically in to. Be- 



cause the superfluid density can be directly related to the 
helicity modulus, it is tempting to describe the strongly 
disordered vortex glass phase by numerically studying 
the gauge glass model, although certain aforementioned 
provisions need to be considered: the gauge glass model 
does not include a field anisotropy and it does not have 
transverse screening. We use our results for the helic- 
ity modulus in the gauge glass to discuss the superfluid 
density in the vortex glass, and thus, test if the super- 
fluid density vanishes in the thermodynamic limit or if 
it stays finite when the system is in equilibrium (ui = 0). 
Our results from equilibrium Monte Carlo simulations of 
the gauge glass for the helicity modulus show that for the 
system sizes studied the superfluid density remains finite 
in the thermodynamic limit. 

The paper is structured as follows: In Sec.|TT]we intro- 
duce the model, observables, and the numerical method 
used. In Sec. IIIII we present our results, followed by con- 
cluding remarks in Sec. IIVI 



II. MODEL, OBSERVABLES, AND 
NUMERICAL DETAILS 

The Hamiltonian of the gauge glass is given by 



w = - E J K 



Aij), 



(1) 



(i>J> 



where the sum ranges over nearest neighbors on a cubic 
lattice in three space dimensions of size N — L 3 , and 
4>i represent the angles of the XY spins. For the gauge 
glass, Jij = J V i, j 'in Eq. (TT]) and in this work we set 
J = 1. Periodic boundary conditions (PBCs) are ap- 
plied. The Aij represents the line integral of the vector 
potential between sites i and j, and we chose the Aij from 
a uniform distribution in [0, 2ir] but with the constraint 
In three space dimensions, the model 



that Ai 
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exhibits a finite-temperature spin-glass transition 1 - 1 - at 
T c « 0.47 ± 0.03. Note that when screening is added 
to the Hamiltonian in Eq. (fT]), the transition to a glass 
phase is lostj£ 

The helicity modulus Y is the second-order derivative 
of the free energy F with respect to an infinitesimal twist 
to the boundaries^ i.e., Y = d 2 F/dG 2 \ e ^ , 

Y = j2 [E< cos ( A *)> - ^E< sin ( A *) sin ( A ^)> 



-(sin(A l ))(sin(A J ))]. (2) 

In Eq. ((2]), (■ • • ) represents a thermal average and Aj = 
4>i — 4>i+x — Au +X , where the twist is performed along the 
x axis.— Because a disorder average of the free energy F 
is independent of the twist angle, we expect that 



[Y] av = (PBC), 



(3) 



where [■ • ■ ] av represents a disorder average. This is not 
necessarily the case for the typical (median) value [V]typ- 
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If we study the distribution P(Y) of Y, then we expect 
a large positive contribution due to spin waves. In ad- 
dition, abrupt changes in the slope of the free energy F 
due to sudden vortex rearrangements in certain disorder 
configurations will generate rare, very negative values of 
the helicity modulus Y, thus generating a long negative 
tail in the distribution^ 6 - (see Fig. EJ. This immediately 
poses the question of whether the different moments of 
the distribution are properly defined. Especially for low 
temperatures T <C T c , we expect the tail of the distribu- 
tion to be most pronounced due to the low probability of 
vortex rearrangements. 

In order to better quantify the effects of the tail of the 
distribution, we empirically fit the data to a generalized 
extreme-value distribution / 6 ! 17 which, in general, is the 
limiting distribution of the maxima of a sequence of in- 
dependent and identically distributed random variables. 
Because of this property, generalized extreme-value dis- 
tributions are used as an approximation to model the 
maxima of long (finite) sequences of random variables 
that can be, for example, found in the analysis of stock 
market data or time series in flat-histogram methods . 28 i 29 
Note that Bertin and Clusel 3 ^ have shown that sums 
of correlated variables also yield limiting distributions, 
which can be described by extreme- value distributions. 
Therefore an underlying extreme process is not necessary 
to obtain an extreme-value limiting distribution. In ad- 
dition to the aforementioned applications, extreme-value 
distributions have recently been used to study different 
problems in other areas of science ranging from physics 
and astronomy to biology applications. In the present 
context, the use of extreme-value distributions is moti- 
vated by the extremely skewed shape of the helicity mod- 
ulus distribution as well as the strong fluctuations found 
in the mean and the standard deviation. 

The cumulative (integrated) generalized extreme- value 
distribution is given by 



Hc(x) = exp 



- 1 + C 



x — fl 



-vr 



(4) 



where /i is the "location parameter" (most probable 
value), (3 the "scale parameter" (standard deviation of 
the peak), and £ is the "shape parameter" which indi- 
cates the asymptotic behavior of the tail of the distribu- 
tion. The generalized extreme-value distribution is then 
given by 



H(x) 



dH c {x) 
dx 



(5) 



where Hq(x) is given by Eq. (gj). In Eq. ©, when £ < 
we obtain a thin-tailed (decay faster than exponential) 
Weibull distribution. If £ = 0, we obtain a Gumbcl 
extreme- value distribution (exponential tail), whereas if 
£ > 0, we obtain a fat-tailed Frechet distribution where 
the tail falls off slower than an exponential. This has the 
effect that when £ > 0, the mth moment of the Frechet 
distribution exists only if £ < 1/m; i.e., if £ > 1/2, the 



variance of the distribution is not defined, and if £ > 1 
the mean is not properly defined either. Here, we expect 
that the distribution of the helicity modulus is well fitted 
by an extreme- value distribution, i.e., 



P(Y)=H(-Y) 



(6) 



where H(x) is given by Eq. ([5]). 

Finally, the superfluid density p s is related to the he- 
licity modulus via 26 



Ps = L^ d Y, 



(7) 



where d represents the space dimension. In the present 
case, d = 3 and thus p s = L~ Y Y . Therefore, by investi- 
gating the properties of the helicity modulus of the gauge 
glass, we can also obtain information about the superfluid 
density of the model. 

The simulations are performed using the parallel tem- 
pering Monte Carlo methodi 31 ' 32 To ensure that the 
data are in equilibrium, we perform a logarithmic bin- 
ning of the moments of different observables (such as the 
current s 11 ' 26 ^ 3 and the helicity modulus defined above) 
and ensure that the data are independent of Monte Carlo 
steps for the last three binsi 3 ^ In addition, we ensure 
that the acceptance probabilities for the parallel temper- 
ing moves are always greater than 0.3. In order to speed 
up the simulations, we discretize the angles of the XY 
spins to N$ = 512 values, a number large enough to 
avoid crossover effects to other models^ Finally, to en- 
sure that the local spin flips are accepted often enough at 
low temperatures (~ 20%), we select the proposed new 
angle for a spin from a temperature-dependent window 
around the original angle, i.e., A(f> cx T with a minimum 
angle for low temperatures of 5°. The simulation pa- 
rameters are presented in Table HI Note that to avoid 
statistical bias in the computation of the helicity modu- 
lus, four replicas of the system have been used for each 
temperature simulated. 

The data are fitted to a generalized extreme- value dis- 
tribution using the statistics package R (see Refs.lT^ and 
36h using the f gev optimization routine in the evd pack- 
age. 



III. RESULTS 

Figure [2] shows the distribution of the helicity modulus 
in the disorder for different system sizes at T = 0.409, 
which is below the known value of T c = 0.47(3) ^ The 
data show a peak for Y > and a long tail that extends 
to large negative values, in agreement with the discussion 
in Sec.|TTJ The behavior of the tail can be studied in more 
detail in a semilogarithmic plot, see Fig. [3J For T < T c , 
the tail decays slower than exponential, thus illustrating 
the possibility of moments of the distribution not being 
properly defined. In contrast, for T > T c , the tail of 
the distribution decays faster than exponential, as can 
be seen, for example, for T = 0.791 in Fig. 2J 
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FIG. 2: (Color online) Distribution of the helicity modulus 
P(Y) for different system sizes L at T = 0.409 < T c . The 
data show a pronounced peak close to Y > (see inset) and 
a long tail for Y < which grows with increasing system size, 
thus suggesting that the data are extreme- value distributed. 
The dashed lines are guides for the eyes. 
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FIG. 3: (Color online) Data for the distribution of the helicity 
modulus P(Y) shown in Fig. [2] in a linear-logarithmic plot. 
The curvature in the tails indicates that the data are Frechet 
distributed with a positive shape parameter £ for T = 0.409 < 
T c . The width of the tail increases slightly with increasing 
system size L. The error bars are not shown for clarity. 



TABLE I: Parameters of the simulation. N sw is the total 
number of Monte Carlo sweeps, T m i n is the lowest temperature 
simulated, and Nt is the number of temperatures used in the 
parallel tempering method (Ref. [37h for each system size L. 
The highest temperature simulated is T max = 0.947 S> T c . 
For each system size, jV sa = 10 4 samples have been simulated 
in order to probe the tails of the distribution in detail. 



L iVsw Tlnin Nt_ 

3 6.0 x 10 3 0.050 53 

4 2.0 x 10 4 0.050 53 

5 6.0 x 10 4 0.050 53 

6 2.0 x 10 5 0.050 53 
8 1.2 x 10 6 0.050 53 
12 1.2 x 10 6 0.303 38 
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FIG. 4: (Color online) Linear-logarithmic plot of the distri- 
bution of the helicity modulus P(Y) for T = 0.791 > T c . The 
tails of the distribution decay faster than exponential, thus 
suggesting that the data are well described by a thin-tailed 
Weibull distribution with £ < 0. The error bars are not shown 
for clarity. 



In Fig. [5j we show the average helicity modulus as a 
function of temperature for different system sizes. Fol- 
lowing previous arguments, we expect [Y] av = 0. This is 
the case for high enough temperatures, yet for low tem- 
peratures the data fluctuate strongly with increasing sys- 
tem size due to the thick tails of the distributions, even 
for a large number of disorder realizations. Because we 
expect [y]av = 0, the data show no evidence of the phase 
transition; i.e., there is no crossing of the data for differ- 
ent system sizes as one would expect from a finite-size 
scaling analysis. This is not the case for the typical value 
(median) of the distribution (Fig. 0, where the data are 
clearly nonzero and cross approximately at the known 
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FIG. 5: (Color online) Average helicity modulus [y] avc as a 
function of temperature for several system sizes. The average 
should be zero within error bars for all temperatures. The 
data show that this is the case for T > T c , whereas for T < 
T c , due to the fat tails of the distributions (see Fig. [3J, the 
fluctuations in the mean as well as the standard deviation are 
extremely large. Note that the data show no signature of the 
phase transition at T c « 0.47 (shaded area). 
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FIG. 6: (Color online) Typical (median) helicity modulus 
[y]ty P as a function of temperature for different system sizes. 
The data are clearly nonzero in equilibrium and display the 
correct finite-size scaling; i.e., the data cross at T c » 0.45 - 
0.50 (shaded area). 



value of the transition temperature. 



In order to improve on these results and properly ac- 
count for the long tails in the distribution of the helicity 
modulus, we empirically fit the data for the distributions 
at different temperatures and system sizes to a general- 
ized extreme- value distribution [Eq. Q]; these are the 
solid lines in Figs. [3] and 2] The data are well repre- 
sented by the fit, although in the tails the fitting routine 
slightly underestimates the tail at extremely low temper- 
atures; however, this does not bias the results because 
the location parameter is mainly determined by the peak 
in the distribution. 

We first study the shape of the helicity modulus dis- 
tributions via the shape parameter £y I n Fig- 13 we 
show the shape parameter as a function of temperature 
for the different system sizes simulated. For T > T c , 
the shape parameter is negative for all L, thus indicat- 
ing that when we are above the transition temperature 
the data are well described via a thin-tailed Weibull dis- 
tribution. This is not the case for T < T c , where the 
shape parameter increases for all system sizes when the 
temperature decreases to zero. In particular, already for 
small system sizes L w 8 and T w 0.10, the shape pa- 
rameter exceeds 1/2; i.e., the variance of the distribution 
is ill-defined. This explains the diverging fluctuations of 
the average helicity modulus in Fig. [5] for T — > 0. By 
extrapolating the data, we expect that for T « 0.10 and 
L > 15, the mean will not be properly defined any longer 
(since, then, £ > 1). 

From a physical standpoint, a thin-tailed distribution 
for T > T c suggests that vortex jumps, while they might 
occur, tend to be exponentially suppressed. This is not 
the case for T < T c , where the tail of the distribution re- 
mains "fat;" i.e., arbitrarily large vortex rearrangements 
can occur. Note that for T — > oo, we expect the distri- 
bution in the thermodynamic limit to become Gaussian, 
which corresponds to £ = —1/2. 

Because in an extreme-value distribution the location 
parameter fi represents the most probable value (similar 
to the median), we show in Fig. [5] the location parameter 
of the helicity modulus distribution, //y, as a function of 
temperature for different system sizes. The data intersect 
cleanly at T = 0.48(2), thus showing clear evidence of 
the glass transition. The fact that the data intersect so 
cleanly suggests that the helicity modulus distribution 
can be well described by extreme- value distributions. In 
the inset of Fig. [5J we show a finite-size scaling plot of 
the data according to 

li Y ~M[LV v (T-T c )], (8) 

where M is a scaling function and v is the critical expo- 
nent of the correlation length^ The best scaling collapse 
is determined by a nonlinear minimization routine 3 -* us- 
ing the software package R. 18 ' 36 The data scale well and 
we estimate T c = 0.48(2) and v — 1.62(20), in agree- 
ment with previous estimates of the critical exponents 
by Olson and Young computed from a finite-size scal- 
ing of the currents [T c = 0.47(3) and v = 1.39(20), see 
Ref. [ll| 4^ Note that a recently introduced extended seal- 
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FIG. 7: (Color online) Shape parameter fy of the helicity 
modulus distribution as a function of temperature for differ- 
ent system sizes. Note that £y < for T > T c (shaded re- 
gion); i.e., the distribution is thin-tailed and the moments of 
the distribution are properly defined. For T < T c , the shape 
parameter is positive; i.e., the distribution is fat tailed. Al- 
ready for T « 0.10 and L > 8, the standard deviation of the 
distribution is not properly defined (£y > 1/2, marked by the 
upper horizontal dashed line). This explains the strong fluc- 
tuations observed in the average helicity presented in Fig. [S] 



FIG. 8: (Color online) Location parameter of the helicity 
modulus distribution /iy as a function of temperature for dif- 
ferent system sizes. The data are nonzero and cross cleanly 
at T c = 0.48(2). The inset shows a finite-size scaling analysis 
of the data according to Eq. ©. L — 3 has been omitted 
in the finite-size scaling analysis due to strong corrections to 
scaling. The observed scaling of the data suggests that gener- 
alized extreme-value distributions describe the helicity mod- 
ulus distribution well. The vertical dashed line in the main 
panel at T — 0.48 is a guide for the eyes. 



ing scheme might yield slightly different estimates for 
the critical exponents; however, we expect the location 
of the transition to remain unchanged. 

Since the average helicity modulus is zero and the lo- 
cation parameter shows the correct critical behavior, we 
compute the superfluid density p° c from Eq. by re- 
placing Y with fiy and plot it as a function of system 
size L for different temperatures in Fig. [9] The data 
decay with increasing system size for all temperatures 
below and above the critical point and are well fitted by 
a phenomenological constant+exponential behavior, i.e., 
p l ° c (L) = a + 6cxp(-cL) (dashed lines in Fig. [5]) with 
fitting probabilities^ Q ~ 0.3 for T -> OM Therefore 
our results suggest that p x ° c {L = oo) is nonzero for the 
system sizes studied; for reference, Fisher et al. predict 
p l ° c ~ a + bj log(w/wo) x where ujq is a characteristic fre- 
quency and x an exponent, and Korshuno\*i^ predicts 
a = 0. We emphasize that our results do not rule out 
the scenario as predicted by Korshunov^ where the su- 
perfluid density is zero in equilibrium: The work of Ko- 
rshunov is based on the existence of a hierarchical dis- 
tribution of metastable states, and it is unclear if such 
an approximation would be valid for the gauge glass. In 
addition, large system sizes, which are numerically inac- 
cessible, are required to probe a logarithmic behavior in 
detail. 



The advantage of fitting the data with extreme-value 
distributions is also illustrated when studying the super- 
fluid density: While the data for the typical estimate also 
display qualitatively the same behavior as the data com- 
puted from the location parameter, the fluctuations are 
considerably larger for the former (not shown) and fitting 
the data to extract any limiting behavior is impossible. 



IV. CONCLUSIONS 

We show that the distribution of the helicity modulus 
in the three-dimensional Gauge glass is well described by 
an extreme- value distribution and demonstrate that the 
typical value of the helicity modulus is nonzero, exhibit- 
ing critical scaling, whereas the behavior of the average 
value is an artifact of the boundary conditions: the aver- 
age estimate is zero and displays strong fluctuations be- 
low the transition temperature. These fluctuations can 
be explained by fitting the data to an extreme- value dis- 
tribution and by studying the shape parameter of the 
distribution. Below T c the data are fat tailed, and thus, 
different moments of the distribution are ill defined. 

In addition, we have computed the superfluid density 
from the location parameter of the extreme-value dis- 
tribution and shown that it is finite in the thermody- 
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FIG. 9: (Color online) Superfluid density p^ 00 computed from 
the location parameter py as a function of system size for 
different temperatures. The dashed lines are fits according to 
pi oc (L) =a + bexp(-cL). 

namic limit. In this work, the equilibrium properties of 
the gauge glass [u> — > limit before L — > oo limit] have 
been studied. Another question of interest is if the lim- 
its L — > oo and to — > are interchangeable. Assuming 



the gauge glass to be a good effective model of the vortex 
glass phase in high-temperature superconductors, our re- 
sults for the gauge glass imply that the superfluid density 
of the vortex glass is finite in the thermodynamic limit 
when the system is in equilibrium. Much larger system 
sizes are required to answer this question beyond reason- 
able doubt. 

Finally, it might be of interest to perform similar stud- 
ies on improved models of vortex glasses, such as the 
recently introduced model from Ref. [H which includes 
the disorder via the couplings Jij between the XY spins 
[see Eq. (JTJ)]. Furthermore, the effects of screening on 
the gauge glass should be revisited at very low temper- 
atures. We hope that this work will find application to 
other studies of observable distributions in physics where 
fat tails are involved. 
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